function [phasor_same, phasor_cross, phasor_perp, phasor_paral] = ...
get_reflection_coeff_layered (elev_incid, perm_bottom, perm_top, perm_middle, sfc, opt)
    if isfieldempty(opt, 'layered_reflection_code')
        if isfield(opt, 'code') && ~isempty(opt.code)
            opt.layered_reflection_code = opt.code;
        else
            opt.layered_reflection_code = 'zavorotny2';
        end
    end
    switch lower(opt.layered_reflection_code)
    case 'orfanidis'
        code = @mymultidiel;
    case {'zavorotny', 'valery'}
        code = @my_refl_mlm;
    case {'zavorotny2', 'valery2'}
        code = @my_refl_mlm2;
    otherwise
        error('snr:get_reflection_coeff_layered:badCode', ...
            'Unkwnown code "%s".', char(opt.layered_reflection_code));
    end
    %code  % DEBUG
    theta = 90 - elev_incid;
    if isempty(perm_bottom),  perm_bottom = perm_middle(end);  end
    perm = [perm_top(:); perm_middle(:); perm_bottom(:)];
    assert(numel(perm) == (numel(sfc.thickness)+2))
    [phasor_perp, phasor_paral] = code(sfc.thickness, perm, opt.frequency, theta);
    [phasor_same, phasor_cross] = get_reflection_coeff_linear2circ (phasor_perp, phasor_paral);
end

%!test
%! % calculate complex relative permittivity:
%! depth = 1e-2*[0; 10; 20; 30];  % interfaces
%! thickness = diff(depth);
%! depth2 = depth(1:end-1,:) + thickness./2;  % layer averages
%! perm = [...
%!   get_permittivity('dry ground');
%!   get_permittivity('medium ground');
%!   get_permittivity('wet ground');
%! ];
%! 
%! %% define auxiliary input:
%! freq = get_gps_const('L2', 'frequency');
%! theta = 0;
%! theta = 10;
%! theta = linspace(0, 90, 10)';
%! elev = 90 - theta;
%! %theta = linspace(0, 90, 10)';
%! 
%! %% test against Fresnel reflection: single interface (two halfspaces)
%! perm_bottom = perm(1)
%! perm_middle = []
%! perm_top = 1
%! [refl_horz_frn, refl_vert_frn] = get_reflection_coeff_linear (elev, perm_bottom, perm_top)
%! [refl_horz_fr2, refl_vert_fr2] = fresnel (sqrt(perm_top), sqrt(perm_bottom), theta, 'std')
%! [ignore, ignore, refl_horz_ewa, refl_vert_ewa] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',[]), struct('code','orfanidis', 'frequency',freq));
%! [ignore, ignore, refl_horz_val, refl_vert_val] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',[]), struct('code','valery', 'frequency',freq));
%! [ignore, ignore, refl_horz_va2, refl_vert_va2] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',[]), struct('code','valery2', 'frequency',freq));
%!   temp = [...
%!    refl_horz_frn ...
%!    refl_horz_fr2 ...
%!    refl_horz_ewa ...
%!    refl_horz_val ...
%!    refl_horz_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_horz_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -100*eps)
%!   temp = [...
%!    refl_vert_frn ...
%!    refl_vert_fr2 ...
%!    refl_vert_ewa ...
%!    refl_vert_val ...
%!    refl_vert_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_vert_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -100*eps)
%! 
%! %% test against Fresnel reflection: single layer (two interfaces, two halfspaces) yet with bottom halfspace similar as the one layer (i.e., with no contrast at the bottom interface):
%! perm_bottom = perm(1)
%! perm_middle = perm(1)
%! perm_top = 1
%! [refl_horz_frn, refl_vert_frn] = get_reflection_coeff_linear (elev, perm_middle, perm_top)
%! [refl_horz_fr2, refl_vert_fr2] = fresnel (sqrt(perm_top), sqrt(perm_middle), theta, 'std')
%! [ignore, ignore, refl_horz_ewa, refl_vert_ewa] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)), struct('code','orfanidis', 'frequency',freq));
%! [ignore, ignore, refl_horz_val, refl_vert_val] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)), struct('code','valery', 'frequency',freq));
%! [ignore, ignore, refl_horz_va2, refl_vert_va2] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)), struct('code','valery2', 'frequency',freq));
%!   temp = [...
%!    refl_horz_frn ...
%!    refl_horz_fr2 ...
%!    refl_horz_ewa ...
%!    refl_horz_val ...
%!    refl_horz_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_horz_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -10*eps)
%!   temp = [...
%!    refl_vert_frn ...
%!    refl_vert_fr2 ...
%!    refl_vert_ewa ...
%!    refl_vert_val ...
%!    refl_vert_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_vert_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -100*eps)
%! 
%! %% test against Fresnel reflection: single thick layer with contrast w.r.t. bottom half-space:
%! %theta = 0;
%! %clear my_refl_mlm refl_mlm
%! %perm_bottom = perm(1)  % WRONG!
%! perm_bottom = perm(2)
%! perm_middle = perm(1)
%! perm_top = 1
%! [refl_horz_frn, refl_vert_frn] = get_reflection_coeff_linear (elev, perm_middle, perm_top)
%! [refl_horz_fr2, refl_vert_fr2] = fresnel (sqrt(perm_top), sqrt(perm_middle), theta, 'std')
%! [ignore, ignore, refl_horz_ewa, refl_vert_ewa] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)*1e6), struct('code','orfanidis', 'frequency',freq));
%! [ignore, ignore, refl_horz_val, refl_vert_val] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)*1e6), struct('code','valery', 'frequency',freq));
%! [ignore, ignore, refl_horz_va2, refl_vert_va2] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!   perm_middle, struct('thickness',thickness(1)*1e6), struct('code','valery2', 'frequency',freq));
%!   assert(~all(isnan(refl_vert_ewa(:))))
%!   assert(~all(isnan(refl_horz_ewa(:))))
%!   assert(~all(isnan(refl_vert_val(:))))
%!   assert(~all(isnan(refl_horz_val(:))))
%!   assert(~all(isnan(refl_vert_va2(:))))
%!   assert(~all(isnan(refl_horz_va2(:))))
%!   temp = [...
%!    refl_horz_frn ...
%!    refl_horz_fr2 ...
%!    refl_horz_ewa ...
%!    refl_horz_val ...
%!    refl_horz_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_horz_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -10*eps)
%!   temp = [...
%!    refl_vert_frn ...
%!    refl_vert_fr2 ...
%!    refl_vert_ewa ...
%!    refl_vert_val ...
%!    refl_vert_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_vert_frn)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -100*eps)
%! 
%! %% test against Fresnel reflection: one layer with increasing thickness:
%! thickness_all = linspace(0,1e3)';
%! %perm_bottom = perm(1)  % WRONG!
%! perm_bottom = perm(2)
%! perm_middle = perm(1)
%! perm_top = 1
%! refl_horz_diff = [];
%! refl_vert_diff = [];
%! for i=1:numel(thickness_all)
%!   %i  % DEBUG
%!   [refl_horz_frn, refl_vert_frn] = get_reflection_coeff_linear (elev, perm_middle, perm_top);
%!   [refl_horz_fr2, refl_vert_fr2] = fresnel (sqrt(perm_top), sqrt(perm_middle), theta, 'std');
%!   [ignore, ignore, refl_horz_ewa, refl_vert_ewa] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness_all(i)), struct('code','orfanidis', 'frequency',freq));
%!   [ignore, ignore, refl_horz_val, refl_vert_val] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness_all(i)), struct('code','valery', 'frequency',freq));
%!   [ignore, ignore, refl_horz_va2, refl_vert_va2] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness_all(i)), struct('code','valery2', 'frequency',freq));
%!   temp = [...
%!    refl_horz_frn ...
%!    refl_horz_fr2 ...
%!    refl_horz_ewa ...
%!    refl_horz_val ...
%!    refl_horz_va2 ...
%!   ];
%!   temp = bsxfun(@minus, temp, refl_horz_frn);
%!   refl_horz_diff = cat(3, refl_horz_diff, temp);
%!   temp = [...
%!    refl_vert_frn ...
%!    refl_vert_fr2 ...
%!    refl_vert_ewa ...
%!    refl_vert_val ...
%!    refl_vert_va2 ...
%!   ];
%!   temp = bsxfun(@minus, temp, refl_vert_frn);
%!   refl_vert_diff = cat(3, refl_vert_diff, temp);
%!   %pause
%! end
%! % temp = permute(real(refl_horz_diff(:,4,:)), [1,3,2]);
%! %   figure, imagesc(elev, thickness_all*100, temp), colorbar
%! refl_horz_diff2 = squeeze(refl_horz_diff(elev==90,5,:));
%! refl_vert_diff2 = squeeze(refl_vert_diff(elev==90,5,:));
%! %refl_horz_diff2 = squeeze(refl_horz_diff(elev==90,4,:));
%! %refl_vert_diff2 = squeeze(refl_vert_diff(elev==90,4,:));
%! %refl_horz_diff2 = squeeze(refl_horz_diff(elev==90,3,:));
%! %refl_vert_diff2 = squeeze(refl_vert_diff(elev==90,3,:));
%!   figure
%!     hold on
%!     plot(real(refl_horz_diff2), thickness_all*100, '.-r', 'LineWidth',1)
%!     plot(real(refl_vert_diff2), thickness_all*100, '+-r', 'LineWidth',1)
%!     plot(imag(refl_horz_diff2), thickness_all*100, '.-b', 'LineWidth',1)
%!     plot(imag(refl_vert_diff2), thickness_all*100, '+-b', 'LineWidth',1)
%!     set(gca, 'YDir','reverse')
%!     grid on
%!     xlabel('Layered minus Fresnel (magnitude of complex difference)')
%!     ylabel('Layer thickness (cm)')
%!     legend({'Real Horizontal','Real Vertical', 'Imag Horizontal','Imag Vertical'}, 'Location','SouthWest');
%!   % nothing to test.
%! 
%! %% general case:
%! perm_bottom = []
%! perm_middle = perm
%! perm_top = 1
%! thickness = [0.1; 0.2; 0.05];
%!   [ignore, ignore, refl_horz_ewa, refl_vert_ewa] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness), struct('code','orfanidis', 'frequency',freq));
%!   [ignore, ignore, refl_horz_val, refl_vert_val] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness), struct('code','valery', 'frequency',freq));
%!   [ignore, ignore, refl_horz_va2, refl_vert_va2] = get_reflection_coeff_layered (elev, perm_bottom, perm_top, ...
%!     perm_middle, struct('thickness',thickness), struct('code','valery2', 'frequency',freq));
%!   temp = [...
%!    refl_horz_ewa ...
%!    refl_horz_val ...
%!    refl_horz_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_horz_ewa)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -10*eps)
%!   temp = [...
%!    refl_vert_ewa ...
%!    refl_vert_val ...
%!    refl_vert_va2 ...
%!   ]
%!   temp = bsxfun(@minus, temp, refl_vert_ewa)
%!   temp(isnan(temp)) = 0;
%!   myassert(temp, 0, -10*eps)


